********************************************************************************
*************************** Do File (5): Table 4 *******************************
********************************************************************************

clear
matrix drop _all
set more off

use "$data/EPH_1015_format.dta", clear

destring occupation occupation_unemp occup_all, replace

*** Keep only observations that have all controls used

keep if (occupation !=. & underemp != . & cont_pens !=. & domwk_all !=. & linc_mjob_base08 !=. & lwagehr_mjob_base08 !=. & linc_total_base08 != . & lhours_mjob !=. & age !=. & hhsize !=. & tenure !=. & msa !=. & educyr !=. ///
& attsch_ever !=. & native !=. & lit !=. & gender !=.) | (occupation_unemp != . & age !=. & hhsize !=. & msa !=. & educyr !=. & attsch_ever !=. & native !=. & lit !=. & gender !=.)

** Keeps only women to run the analysis using observations from only one household member
keep if gender == 0

********************************************************************************
********************************************************************************

*** Locals for controls

local occup occupation
local base_controls "msa year"
local controls "age age2 hhsize lit native attsch_ever educyr educyr2 i.marstat i.dec_pcfaminc"

*** Label the variables as they appear in the table

label var mean_active "Labor force participation"
label var share_formal "Share registered"
label var lhhd_hours_mjob "Hours of work per week"
label var lhhd_labinc_base08 "Labor income per month"

********************************************************************************
********************************************************************************
********************************************************************************

* Run regression for pension contribution to restrict the sample to use

qui reghdfe cont_pens domwk_all treat_dwall `controls' if ${ctrl_group} == 1 & unemployed == 0, absorb(`base_controls' `occup') vce(cluster msa)
gen sample_reg = 1 if e(sample)


********************************************************************************
********************************************************************************

** Households with spouse or children

sort time CODUSU nro_hogar pid

duplicates drop time CODUSU nro_hogar, force


* First run all the regressions to get the p-values and obtain adjusted p-values

foreach var in mean_active share_formal lhhd_hours_mjob lhhd_labinc_base08 {
	
	qui reghdfe `var' hhd_domwk treat_hhddw if sample_reg == 1 & (haschild1625 == 1 | marstat == 2), absorb(`base_controls') vce(cluster msa)

	capture matrix list uP_treat
	
	if _rc != 0 {
	    
		matrix define uP_treat = 2*ttail(e(N),abs(_b[treat_hhddw]/_se[treat_hhddw]))
		
	}
	
	else {
	    
		matrix uP_treat = uP_treat, 2*ttail(e(N),abs(_b[treat_hhddw]/_se[treat_hhddw]))
	}
	
}

* Generate the adjusted p-values and put them in a matrix that will be exported with the regression estimates

matrix uP_treat = uP_treat'
svmat uP_treat, names(unad_p)

qqvalue unad_p1, method(hochberg) qvalue(hochbergP)

qui sum unad_p1
mkmat unad_p1 hochbergP if _n <= r(N), mat(adjPval)
matrix adjPval = adjPval[1..r(N),2]

** Rerun all the regressions, including the adjusted p-values

local i = 1

local varlabel: variable label mean_active
qui reghdfe mean_active hhd_domwk treat_hhddw if sample_reg == 1 & (haschild1625 == 1 | marstat == 2), absorb(`base_controls') vce(cluster msa)
sum mean_active if hhd_domwk == 1 & treat == 0 & e(sample)
local meanvar = r(mean)
local qval = round(adjPval[`i',1], 0.001)
outreg2 using "$tables/Table_4", replace excel keep(treat_hhddw) nocons dec(3) label ctitle(`varlabel') addstat(Mean dependent variable,`meanvar', q-value, `qval') ///
addtext(Controls, No, Year Fixed Effects, Yes, Occupation Fixed Effects, No, MA Fixed Effects, Yes, Number of clusters, "`e(N_clust)'")

local i = `i' + 1


local varlabel: variable label share_formal
qui reghdfe share_formal hhd_domwk treat_hhddw if sample_reg == 1 & (haschild1625 == 1 | marstat == 2), absorb(`base_controls') vce(cluster msa)
sum share_formal if hhd_domwk == 1 & treat == 0 & e(sample)
local meanvar = r(mean)
local qval = round(adjPval[`i',1], 0.001)
outreg2 using "$tables/Table_4", append excel keep(treat_hhddw) nocons dec(3) label ctitle(`varlabel') addstat(Mean dependent variable,`meanvar', q-value, `qval') ///
addtext(Controls, No, Year Fixed Effects, Yes, Occupation Fixed Effects, No, MA Fixed Effects, Yes, Number of clusters, "`e(N_clust)'")

local i = `i' + 1

foreach var in hhd_hours_mjob hhd_labinc_base08 {
	
	local varlabel: variable label l`var'
	qui reghdfe l`var' hhd_domwk treat_hhddw if sample_reg == 1 & (haschild1625 == 1 | marstat == 2), absorb(`base_controls') vce(cluster msa)
	sum `var' if hhd_domwk == 1 & treat == 0 & e(sample)
	local meanvar = r(mean)
	local qval = round(adjPval[`i',1], 0.001)
	outreg2 using "$tables/Table_4", append excel keep(treat_hhddw) nocons dec(3) label ctitle(`varlabel') addstat(Mean dependent variable,`meanvar', q-value, `qval') ///
	addtext(Controls, No, Year Fixed Effects, Yes, Occupation Fixed Effects, No, MA Fixed Effects, Yes, Number of clusters, "`e(N_clust)'")
	
	local i = `i' + 1
	
}